//global outpath "C:/Dropbox/projects/groundwater/replication"
global outpath E:/Work/Dropbox/projects/groundwater/replication

cap cd "${outpath}"









/************************************************************************
*************************************************************************
*************************************************************************

Main Text

*************************************************************************
*************************************************************************
************************************************************************/



/************************************************************************
*************************************************************************
Figure 1
*************************************************************************
************************************************************************/



/************************************************************************
Water Levels
************************************************************************/

use "dta/crosscut_station_long.dta", clear


foreach var of varlist idw_pmr* {
	gen p50`var' = `var'
	gen p75`var' = `var'
	gen p25`var' = `var'
}

collapse (p25) p25* (p50) p50* (p75) p75*

gen i = 1
reshape long p25idw_pmr p50idw_pmr p75idw_pmr  , i(i) j(yr) string

destring yr, replace
replace yr = 1900 + yr if yr > 20
replace yr = yr + 2000 if yr < 20

gen textlab_med = "Median" if yr == 2002
gen textlab_iqr = "Inter-Quartile Range" if yr == 2010

gen arrow_x1 = 1995
gen arrow_y1 = 4
gen arrow_x2 = 1995
gen arrow_y2 = 7.75

gen arrow_lab = "Greater Water Stress"

#d;
twoway
	(scatter p50idw_pmr yr, connect(direct) lp(solid) lc(black) mc(black) sort(yr))
	(rarea p75idw_pmr p25idw_pmr yr, color(gs8%25) sort lw(none) )
	(scatter p50idw_pmr yr, ms(i) mlab(textlab_med) mlabcolor(black) mlabpos(3) mlabsize(medlarge) )
	(scatter p25idw_pmr yr, ms(i) mlab(textlab_iqr) mlabcolor(gs8) mlabpos(12) mlabsize(medlarge) )
	(pcarrow arrow_y1 arrow_x1 arrow_y2 arrow_x2, lc(navy) mc(navy) mlabel(arrow_lab) mlabcolor(navy) mlabpos(1) mlabangle(-90) mlabsize(medlarge) )
	, 
	
	legend(off) graphregion(color(white))
	xtitle("Year") ytitle("Drill Depth (Meters From Surface, Inverted Scale)") yscale(reverse)
	yline(8, lc(navy) lp(shortdash) ) text(8 2005 "Max Depth of Standard Pumps", place(6) color(navy) )
	;
#d cr

cap mkdir "${outpath}/results"
cap mkdir "${outpath}/results/motivation"
graph export "${outpath}/results/motivation/declining_water.pdf", replace





/************************************************************************
Correlations
************************************************************************/

use "dta/nregs_gp.dta", clear


preserve
	keep if y > 13
	keep if y==14
	gen floordepth = floor(min_nearestdepth) + .5
	collapse rabi_TOTmeandaysworked , by(floordepth)

	#d;
		twoway
		( scatter rabi_TOTmeandaysworked floordepth if floordepth  < 12)
		( lfit rabi_TOTmeandaysworked floordepth if floordepth  < 12) 
		, 
		xline(8, lc(gs8) lp(dash))
		graphregion(color(white))
		xtitle("Drilling Depth (Dry Season, 2014)", size(large) ) ytitle("Average NREGA Days in 2014", size(large) )
		name(g1)
		legend(off)
		;
	#d cr
restore



use "dta/crosscut_station.dta", clear

preserve
	gen floordepth = floor(min_nearestdepth) + .5
	collapse frelAITC , by(floordepth)

	#d;
		twoway
		( scatter frelAITC floordepth if floordepth  < 12)
		( lfit frelAITC floordepth if floordepth  < 12)
		, 
		xline(8, lc(gs8) lp(dash))
		graphregion(color(white))
		xtitle("Drilling Depth (Dry Season, 2014)", size(large) ) ytitle("Average Station-Level Vote Lean Towards Ruling Party", size(large) )
		name(g2)
		legend(off)
		;
	#d cr
restore

graph combine g1 g2, xsize(8) rows(1) graphregion(color(white) fcolor(white))
graph drop g1 g2

cap mkdir "results/suggestive"
graph export "results/suggestive/strawman.pdf", replace





/************************************************************************
*************************************************************************
Figure 2
*************************************************************************
************************************************************************/
/*
See ArcGIS Project File
*/




/************************************************************************
*************************************************************************
Figure 3
*************************************************************************
************************************************************************/


/************************************************************************
NREGA
************************************************************************/


use "dta/nregs_gp.dta", clear
keep if y > 13
rdrobust rabi_TOTmeandaysworked taxi , vce(cluster pid)
gen h = e(h_l)

local aggvar meanidw_pmrb14
sum `aggvar', detail
gen low_water = `aggvar' > r(p50) if `aggvar' < . 
gen high_water = `aggvar' <= r(p50) if `aggvar' < . 
gen overall = 1


reghdfe rabi_TOTmeandaysworked, absorb(i.did i.pcnumber) res(rabi_TOTmeandaysworked_res)


foreach sample of varlist low_water high_water {
	cap graph drop `sample'
	
	preserve
		#d;
			rdplot rabi_TOTmeandaysworked_res taxi if abs(taxi) < h & `sample' == 1, 
			kernel(triangular) p(1) nbins(15) hide genvars 
			;
		#d cr
		
		keep rdplot_*
		duplicates drop
		drop if mi(rdplot_N)
		gen sample = "`sample'"
		
		tempfile `sample'
		save ``sample''
	restore
		
}

clear


append using `low_water' `high_water'


egen tag = tag(sample rdplot_id)



local title_low_water "Water-Stressed"
local title_high_water "Non-Stressed"


sum rdplot_mean_y
local miny = r(min)
local maxy = r(max)

local yscmin = `miny' - mod(`miny',.5)
local yscmax = `maxy' - mod(`maxy',.5) + .5


foreach sample in low_water high_water {
	#d;
	twoway 
		(scatter rdplot_mean_y rdplot_mean_bin if sample=="`sample'" & tag, mc(gs8%25) mlw(none) )
		(lfit rdplot_hat_y rdplot_mean_bin if sample=="`sample'" & rdplot_mean_bin < 0, sort(rdplot_mean_bin) lc(black) )
		(lfit rdplot_hat_y rdplot_mean_bin if sample=="`sample'" & rdplot_mean_bin >= 0, sort(rdplot_mean_bin) lc(black) )
		,
		graphregion(color(white)) legend(off) xtitle("Distance to Ruling Party Majority") name(`sample') 
		xline(0, lp(dash) lc(gs8)) title("`title_`sample''") `ysc`sample''
		ylabel(`yscmin' 0 `yscmax') ymtick(`yscmin'(.5)`yscmax')  ytitle("Dry Season NREGA Days Per Household (Residualized)")
	;
	#d cr
}

graph combine low_water high_water, rows(1) graphregion(color(white) fcolor(white)) ysize(3.5) 
graph drop low_water high_water


cap mkdir "${outpath}/results/nrega"
graph export "${outpath}/results/nrega/nrega_rd_graph.pdf", replace





/************************************************************************
Votes
************************************************************************/

use "dta/crosscut_station.dta", clear
rdrobust frelAITC taxi , vce(cluster pid)
gen h = e(h_l)

/*
Define dummy for whether the water depth---how far you have to dig to get to
water---is above median
*/
sum idw_pmrb14, detail
gen low_water = idw_pmrb14 > r(p50) if idw_pmrb14 < . 
gen high_water = idw_pmrb14 <= r(p50) if idw_pmrb14 < . 
gen overall = 1

reghdfe frelAITC, absorb(i.did i.pcnumber) res(frelAITC_res)


foreach sample of varlist low_water high_water {
	cap graph drop `sample'
	
	preserve
		#d;
			rdplot frelAITC_res taxi if abs(taxi) < h & `sample' == 1, 
			kernel(triangular) p(1)  hide genvars nbins(15)
			;
		#d cr
		
		keep rdplot_*
		duplicates drop
		drop if mi(rdplot_N)
		gen sample = "`sample'"
		
		tempfile `sample'
		save ``sample''
	restore
		
}

clear


append using `low_water' `high_water'


egen tag = tag(sample rdplot_id)



local title_low_water "Water-Stressed"
local title_high_water "Non-Stressed"


sum rdplot_mean_y
local miny = r(min)
local maxy = r(max)

local yscmin = `miny' - mod(`miny',.05)
local yscmax = `maxy' - mod(`maxy',.05) + .05


foreach sample in low_water high_water {
	#d;
	twoway 
		(scatter rdplot_mean_y rdplot_mean_bin if sample=="`sample'" & tag, mc(gs8%25) mlw(none) )
		(lfit rdplot_hat_y rdplot_mean_bin if sample=="`sample'" & rdplot_mean_bin < 0, sort(rdplot_mean_bin) lc(black) )
		(lfit rdplot_hat_y rdplot_mean_bin if sample=="`sample'" & rdplot_mean_bin >= 0, sort(rdplot_mean_bin) lc(black) )
		,
		graphregion(color(white)) legend(off) xtitle("Distance to Ruling Party Majority") name(`sample') 
		xline(0, lp(dash) lc(gs8)) title("`title_`sample''") `ysc`sample''
		ylabel(`yscmin' 0 `yscmax') ymtick(`yscmin'(.05)`yscmax')  ytitle("Ruling Party Vote Share (Residualized)")
	;
	#d cr
}

graph combine low_water high_water, rows(1) graphregion(color(white) fcolor(white)) ysize(3.5)
graph drop low_water high_water

cap mkdir "${outpath}/results/votes"
cap mkdir "${outpath}/results/votes/stationlevel"
graph export "${outpath}/results/votes/stationlevel/rd_graph.pdf", replace




/************************************************************************
*************************************************************************
Figure 4
*************************************************************************
************************************************************************/


/************************************************************************
decollinear
Takes a varlist and returns a list with a non-collinear subset
************************************************************************/
cap program drop decollinear
program decollinear, rclass
	syntax varlist [if]
	
	_rmcoll `varlist' `if'
	local decoll
	foreach term in `r(varlist)' {
		if substr("`term'",1,2) != "o." {
			local decoll `decoll' `term'
		}
	}
	
	return clear
	return local varlist `decoll'
end






/************************************************************************
Vote Shares
************************************************************************/
use "dta/crosscut_station.dta", clear


gen high_water = min_nearestdepth < 8 if !mi(min_nearestdepth)
gen low_water = min_nearestdepth >= 8 if !mi(min_nearestdepth)

keep if min_nearestdist < 5


cap mkdir "estimates"
cap postclose est
postfile est str20 sample bandwidth b cl cu rcl rcu n using "estimates/8mcutoff_votes", replace
	
quietly {
	
	
	foreach p of numlist 8(-.25)1 {
		
			
			
		foreach sample in low_water high_water {
			preserve
				keep if `sample'==1 & abs(min_nearestdepth - 8) < `p'
				rdrobust frelAITC taxi , vce(cluster pid)
				local h = e(h_l)
				
				decollinear INDdist2-INDdist17 INDpcn2-INDpcn37 if abs(taxi) < `h'
				
				rdrobust frelAITC taxi , vce(cluster pid) h(`h') covs(`r(varlist)')
				local b		= _b[RD_Estimate]
				local cu	=  e(ci_r_cl)
				local cl	=  e(ci_l_cl)
				
				local rcu	=  e(ci_r_rb)
				local rcl	=  e(ci_l_rb)
				
				egen tag = tag(pid) if abs(taxi) < e(h_l) & abs(min_nearestdepth - 8) < `p'
				count if tag == 1
				local nclust = r(N)
				drop tag
				
				post est ("`sample'") (`p') (`b') (`cl') (`cu') (`rcl') (`rcu') (`nclust')
			restore
		}
		
		
	}
	
}

postclose est





/************************************************************************
NREGA
************************************************************************/
use "dta/nregs_gp.dta", clear
drop if y==13

gen high_water = min_nearestdepth < 8 if !mi(min_nearestdepth)
gen low_water = min_nearestdepth >= 8 if !mi(min_nearestdepth)

keep if min_nearestdist < 5


cap mkdir "estimates"
cap postclose est
postfile est str20 sample bandwidth b cl cu rcl rcu n using "estimates/8mcutoff_nrega", replace

quietly {
	
	
	foreach p of numlist 8(-.25)1 {
		
			
			
		foreach sample in low_water high_water {
			preserve
				keep if `sample'==1 & abs(min_nearestdepth - 8) < `p'
				rdrobust rabi_TOTmeandaysworked taxi , vce(cluster pid)
				local h = e(h_l)
				
				decollinear INDdist2-INDdist17 INDpcn2-INDpcn37 if abs(taxi) < `h'
				
				rdrobust rabi_TOTmeandaysworked taxi , vce(cluster pid) h(`h') covs(`r(varlist)')
				local b		= _b[RD_Estimate]
				local cu	=  e(ci_r_cl)
				local cl	=  e(ci_l_cl)
				
				local rcu	=  e(ci_r_rb)
				local rcl	=  e(ci_l_rb)
				
				egen tag = tag(pid) if abs(taxi) < e(h_l) & abs(min_nearestdepth - 8) < `p'
				count if tag == 1
				local nclust = r(N)
				drop tag
				
				post est ("`sample'") (`p') (`b') (`cl') (`cu') (`rcl') (`rcu') (`nclust')
			restore
		}
		
		
	}
	
}

postclose est







/************************************************************************
Graphs
************************************************************************/
global intextnotesize medlarge
global rdestcol			maroon
global barheightcol		ltblue

/************
Vote Shares
*************/
use "estimates/8mcutoff_votes", clear

gen rdest_x1 = 3
gen rdest_y1 = .015 
gen rdest_x2 = 1.5
gen rdest_y2 = .015 
gen rdest_lab = "RD Estimate and 95% CI"

gen pointhist_x1 = 5.5
gen pointhist_y1 = 600 
gen pointhist_x2 = 8
gen pointhist_y2 = 600 
gen pointhist_lab = "Bar Height: Clusters in Optimal BW"

local title_low_water "a) Drill Depth Greater than 8m (Stressed)"
local title_high_water "b) Drill Depth Less than 8m (Unstressed)"

local x_low_water 3
local x_high_water 3
foreach sample in low_water high_water {
	
	
	if "`sample'" == "high_water" {
		local xtitles xtitle("Restricting to wells within __m of the 8m Cutoff", size(large) ) xlabel(1(1)8) xscale( range(1 8))
		local arrowlab (pcarrow pointhist_y1 pointhist_x1 pointhist_y2 pointhist_x2, yaxis(2) lc(ltblue%75) mc(ltblue%75) mlabel(pointhist_lab) mlabcolor(${barheightcol}) mlabpos(12) mlabsize(${intextnotesize}) )
	} 
	else {
		local xtitles xtitle("")  xscale( range(1 8) lw(none) )    xlabel(none) 
		local arrowlab (pcarrow rdest_y1 rdest_x1 rdest_y2 rdest_x2, yaxis(1) lc(maroon%50) mc(maroon%50) mlabel(rdest_lab) mlabcolor(${rdestcol}) mlabpos(6)  mlabsize(${intextnotesize}) )
	}
	
	cap graph drop `sample'
	
	#d;
	twoway 
		/*(rspike rcu rcl bandwidth if sample == "`sample'", lc(gs14) lw(thick)) */
		(rspike cu cl bandwidth if sample == "`sample'", lc(maroon%50) ) 
		(scatter b bandwidth if sample == "`sample'", mc(maroon) )
		(bar n bandwidth if sample == "`sample'", color(ltblue%20) barw(.25) yaxis(2) )
		`arrowlab'
		, 
		yline(0, lc(black) lw(thick) ) ylabel(-0.04(.04)0.08, nogrid labcolor(maroon)) ylabel(0(300)900, axis(2) grid labcolor(ltblue) )
		graphregion(color(white))   legend(off) `xtitles' ytitle("")
		 name(`sample')  ytitle("", axis(2) )
		 text(.08 `x_`sample'' "`title_`sample''", color(gs8) placement(3) size(medlarge ) )
		 /*ytitle("Number of Clusters in Optimal BW", axis(2))*/
		/* title("`title_`sample''") ytitle("RD Estimate", size(large) )*/
	;
	#d cr
}

cap graph drop votes
graph combine low_water high_water, rows(2) graphregion(color(white) fcolor(white)) imargin(0 0 0 0) name(votes) title("AITC Vote Lean")
graph drop low_water high_water





/************
Vote Shares
*************/
use "estimates/8mcutoff_nrega", clear

gen rdest_x1 = 2.5
gen rdest_y1 = 4.25
gen rdest_x2 = 1
gen rdest_y2 = 4.25
gen rdest_lab = "RD Estimate and 95% CI"

gen pointhist_x1 = 5.5
gen pointhist_y1 = 600 
gen pointhist_x2 = 8
gen pointhist_y2 = 600 
gen pointhist_lab = "Bar Height: Clusters in Optimal BW"

local title_low_water "a) Drill Depth Greater than 8m (Stressed)"
local title_high_water "b) Drill Depth Less than 8m (Unstressed)"

local x_low_water 3
local x_high_water 3
foreach sample in low_water high_water {
	
	
	if "`sample'" == "high_water" {
		local xtitles xtitle("Restricting to wells within __m of the 8m Cutoff", size(large) ) xlabel(1(1)8) xscale( range(1 8))
		local arrowlab (pcarrow pointhist_y1 pointhist_x1 pointhist_y2 pointhist_x2, yaxis(2) lc(ltblue%75) mc(ltblue%75) mlabel(pointhist_lab) mlabcolor(${barheightcol}) mlabpos(12)  mlabsize(${intextnotesize}) )
	} 
	else {
		local xtitles xtitle("")  xscale( range(1 8) lw(none) )    xlabel(none) 
		local arrowlab (pcarrow rdest_y1 rdest_x1 rdest_y2 rdest_x2, yaxis(1) lc(maroon%50) mc(maroon%50) mlabel(rdest_lab) mlabcolor(${rdestcol}) mlabpos(12)  mlabsize(${intextnotesize}) )
	}
	
	cap graph drop `sample'
	
	#d;
	twoway 
		/*(rspike rcu rcl bandwidth if sample == "`sample'", lc(gs14) lw(thick)) */
		(rspike cu cl bandwidth if sample == "`sample'", lc(maroon%50) ) 
		(scatter b bandwidth if sample == "`sample'", mc(maroon) )
		(bar n bandwidth if sample == "`sample'", color(ltblue%20) barw(.25) yaxis(2) )
		`arrowlab'
		, 
		yline(0, lc(black) lw(thick) ) ylabel(-2(2)6, nogrid labcolor(maroon)) ylabel(0(300)900, axis(2) grid labcolor(ltblue) )
		graphregion(color(white))   legend(off) `xtitles' ytitle("")
		 name(`sample')  ytitle("", axis(2) )
		 text(6 `x_`sample'' "`title_`sample''", color(gs8) placement(3) size(medlarge ) )
		 /*ytitle("Number of Clusters in Optimal BW", axis(2))*/
		/* title("`title_`sample''") ytitle("RD Estimate", size(large) )*/
	;
	#d cr
}

cap graph drop nrega
graph combine low_water high_water, rows(2) graphregion(color(white) fcolor(white)) imargin(0 0 0 0) name(nrega) title("NREGA Days, Dry Season") //xcommon
graph drop low_water high_water




/************
Combine
*************/

graph combine nrega votes , rows(1) graphregion(color(white) fcolor(white)) xsize(7) //imargin(0 0 0 0) //xcommon
graph drop votes nrega


cap mkdir "${outpath}/results/checks"
cap mkdir "${outpath}/results/checks/8m"
graph export "${outpath}/results/checks/8m/cutoff_range_figure.pdf", replace











/************************************************************************
*************************************************************************
Table 1
*************************************************************************
************************************************************************/




/************************************************************************
decollinear
Takes a varlist and returns a list with a non-collinear subset
************************************************************************/
cap program drop decollinear
program decollinear, rclass
	syntax varlist [if]
	
	_rmcoll `varlist' `if'
	local decoll
	foreach term in `r(varlist)' {
		if substr("`term'",1,2) != "o." {
			local decoll `decoll' `term'
		}
	}
	
	return clear
	return local varlist `decoll'
end


/************************************************************************
addnclust
Takes a variable name and calculates
************************************************************************/
cap program drop addnclust
program addnclust
	args var rv h ifstate
	if "`ifstate'" != "" {
		gen SAMPDEF = `ifstate'
	}
	else {
		gen SAMPDEF = 1
	}
	
	egen nmiss = rowmiss(`var' pid `rv')
	egen tag = tag(pid) if nmiss==0 & abs(`rv') < `h' & SAMPDEF==1
	count if tag == 1
	estadd scalar N_clust = r(N)
	drop nmiss tag SAMPDEF
end


/************************************************************************
addconmean
Takes an outcome, a variable name, and a bandwidth to calci;ate amd add
the ``control mean''
************************************************************************/
cap program drop addconmean
program addconmean
	args var rv h ifstate
	if "`ifstate'" != "" {
		gen SAMPDEF = `ifstate'
	}
	else {
		gen SAMPDEF = 1
	}
	
	egen nmiss = rowmiss(`var' pid `rv')
	sum `var' if nmiss==0 & `rv' < 0 & abs(`rv') < `h'/10 & SAMPDEF==1
	estadd scalar conmean = r(mean)
	drop nmiss SAMPDEF
end



/************************************************************************
Votes
************************************************************************/

cap mkdir "${outpath}/results"
cap mkdir "${outpath}/results/votes"
cap mkdir "${outpath}/results/votes/stationlevel"


use "dta/crosscut_station.dta", clear

/*
Define dummy for whether the water depth---how far you have to dig to get to
water---is above median
*/
sum idw_pmrb14, detail
gen low_water = idw_pmrb14 > r(p50) if idw_pmrb14 < . 
gen high_water = idw_pmrb14 <= r(p50) if idw_pmrb14 < . 
gen overall = 1

eststo clear

local out frelAITC

foreach sample of varlist overall low_water high_water {
	//Basic specification
	rdrobust frelAITC taxi if `sample'==1, vce(cluster pid)
	local h = e(h_l)
	estadd local metric = "1-Norm"
	estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
	addnclust `out' taxi `h' "`sample'==1"
	addconmean `out' taxi `h' "`sample'==1"
	eststo
	
	
	//Controlling for fixed-effects (use bandwidth of earlier specification)
	decollinear INDdist2-INDdist17 INDpcn2-INDpcn37 if `sample'==1  & abs(taxi) < `h'

	rdrobust frelAITC taxi if `sample'==1, h(`h') covs(`r(varlist)') vce(cluster pid)
	estadd local distfe = "X"
	estadd local pcfe = "X"
	estadd local metric = "1-Norm"
	estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
	addnclust `out' taxi `h' "`sample'==1"
	addconmean `out' taxi `h' "`sample'==1"
	eststo
}

#d;
esttab , 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) 
sca("N_eff Obs in BW" "N_clust Clusters in BW" "h_l Bandwidth" "pv_rb Robust p-val"   "metric Metric" "distfe District FEs" "pcfe Constituency FEs") 
sfmt("%9.0g" "%9.0g" "%9.3f") replace nonote
mtitles("Overall" "Water-Stressed" "Non-Stressed" "Overall" "Water-Stressed" "Non-Stressed" )
varlabel(RD_Estimate "RD Estimate")
;
#d cr

#d;
esttab using "${outpath}/results/votes/stationlevel/mainresults.tex", 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) 
sca("N Total Obs" "N_eff Obs in BW" "N_clust Clusters in BW" "h_l Bandwidth" "pv_rb Robust p-val"  "distfe District FEs" "pcfe Constituency FEs") 
sfmt("%9.0g" "%9.0g" "%9.0g" "%9.3f") replace nonote
nomtitles noobs
	mgroups(
		"Full Sample" "Water-Stressed" "Non-Stressed",
		pattern(1 0 1 0 1 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
varlabel(RD_Estimate "RD Estimate")
;
#d cr




/************************************************************************
NREGA
************************************************************************/

cap mkdir "results"
cap mkdir "results/nrega"

/************************************************************************
addnclust
Takes a variable name and calculates
************************************************************************/
cap program drop addnclust
program addnclust
	args var rv h ifstate
	if "`ifstate'" != "" {
		gen SAMPDEF = `ifstate'
	}
	else {
		gen SAMPDEF = 1
	}
	
	egen nmiss = rowmiss(`var' pid `rv')
	egen tag = tag(pid) if nmiss==0 & abs(`rv') < `h' & SAMPDEF==1
	count if tag == 1
	estadd scalar N_clust = r(N)
	drop nmiss tag SAMPDEF
end


/************************************************************************
addconmean
Takes an outcome, a variable name, and a bandwidth to calci;ate amd add
the ``control mean''
************************************************************************/
cap program drop addconmean
program addconmean
	args var rv h ifstate
	if "`ifstate'" != "" {
		gen SAMPDEF = `ifstate'
	}
	else {
		gen SAMPDEF = 1
	}
	
	egen nmiss = rowmiss(`var' pid `rv')
	sum `var' if nmiss==0 & `rv' < 0 & abs(`rv') < `h'/10 & SAMPDEF==1
	estadd scalar conmean = r(mean)
	drop nmiss SAMPDEF
end



/************************************************************************
Pooled
************************************************************************/

use "dta/nregs_gp.dta", clear

keep if y > 13
local aggvar meanidw_pmrb14
sum `aggvar', detail
gen low_water = `aggvar' > r(p50) if `aggvar' < . 
gen high_water = `aggvar' <= r(p50) if `aggvar' < . 
gen overall = 1

eststo clear
local out rabi_TOTmeandaysworked
foreach sample of varlist overall low_water high_water {
	quietly {
		
		preserve
			keep if `sample'==1
			

			//Basic specification [we'll use the bandwidth for this in all specifications]
			rdrobust `out' taxi, vce(cluster pid)
			estadd local metric = "1-Norm"
			estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
			addnclust `out' taxi `=e(h_l)'
			addconmean `out' taxi `=e(h_l)'
			eststo
			local modelnames `modelnames' Basic
			local h = e(h_l)


			//Both [Preferred spec]
			decollinear INDdist2-INDdist17 INDpcn2-INDpcn37 if abs(taxi) < `h'
			rdrobust `out' taxi, vce(cluster pid) h(`h') covs(`r(varlist)')
			estadd local distfe = "X"
			estadd local pcfe = "X"
			estadd local metric = "1-Norm"
			estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
			addnclust `out' taxi `h'
			addconmean `out' taxi `h'
			eststo
		restore

	}
}
#d;
esttab, 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) 
sca("N_eff Obs in BW" "N_clust Clusters in BW" "conmean Control Mean" "h_l Bandwidth" "pv_rb Robust p-val"  "metric Metric" "distfe District FEs" "pcfe Constituency FEs") 
sfmt("%9.0g" "%9.0g" "%9.2f" "%9.3f") mtitles("Overall" "Overall" "Water-Stressed" "Water-Stressed" "Non-Stressed"  "Non-Stressed" ) noobs
;
#d cr



#d;
esttab using "${outpath}/results/nrega/pooled.tex", 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) replace nonote
sca("N Total Obs" "N_eff Obs in BW" "N_clust Clusters in BW" "conmean Control Mean" "h_l Bandwidth" "pv_rb Robust p-val"  "distfe District FEs" "pcfe Constituency FEs") 
sfmt("%9.0g" "%9.0g" "%9.0g" "%9.2f" "%9.3f") noobs
nomtitles
	mgroups(
		"Full Sample" "Water-Stressed" "Non-Stressed",
		pattern(1 0 1 0 1 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
varlabel(RD_Estimate "RD Estimate")
;
#d cr






/************************************************************************
*************************************************************************
Table 2
*************************************************************************
************************************************************************/


/************************************************************************
decollinear
Takes a varlist and returns a list with a non-collinear subset
************************************************************************/
cap program drop decollinear
program decollinear, rclass
	syntax varlist [if]
	
	_rmcoll `varlist' `if'
	local decoll
	foreach term in `r(varlist)' {
		if substr("`term'",1,2) != "o." {
			local decoll `decoll' `term'
		}
	}
	
	return clear
	return local varlist `decoll'
end


/************************************************************************
addnclust
Takes a variable name and calculates
************************************************************************/
cap program drop addnclust
program addnclust
	args var rv h ifstate
	if "`ifstate'" != "" {
		gen SAMPDEF = `ifstate'
	}
	else {
		gen SAMPDEF = 1
	}
	
	egen nmiss = rowmiss(`var' pid `rv')
	egen tag = tag(pid) if nmiss==0 & abs(`rv') < `h' & SAMPDEF==1
	count if tag == 1
	estadd scalar N_clust = r(N)
	drop nmiss tag SAMPDEF
end


/************************************************************************
addconmean
Takes an outcome, a variable name, and a bandwidth to calci;ate amd add
the ``control mean''
************************************************************************/
cap program drop addconmean
program addconmean
	args var rv h ifstate
	if "`ifstate'" != "" {
		gen SAMPDEF = `ifstate'
	}
	else {
		gen SAMPDEF = 1
	}
	
	egen nmiss = rowmiss(`var' pid `rv')
	sum `var' if nmiss==0 & `rv' < 0 & abs(`rv') < `h'/10 & SAMPDEF==1
	estadd scalar conmean = r(mean)
	drop nmiss SAMPDEF
end





cap mkdir "${outpath}/results"
cap mkdir "${outpath}/results/checks"
cap mkdir "${outpath}/results/checks/8m"


eststo clear

/************************************************************************
NREGS
************************************************************************/

use "dta/nregs_gp.dta", clear
drop if y==13
/*
Define dummy for whether the water depth---how far you have to dig to get to
water---is above median
*/
gen high_water = min_nearestdepth < 8 if !mi(min_nearestdepth)
gen low_water = min_nearestdepth >= 8 if !mi(min_nearestdepth)

//Stick to 3m around the 8m depth cutoff
keep if abs(min_nearestdepth - 8) < 3

//Keep only polling stations within 5km of a monitoring well
keep if min_nearestdist < 5



local out rabi_TOTmeandaysworked
foreach sample of varlist low_water high_water {
	quietly {
		
		preserve
			keep if `sample'==1
			

			//Basic specification [we'll use the bandwidth for this in all specifications]
			rdrobust `out' taxi, vce(cluster pid)
			estadd local metric = "1-Norm"
			estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
			addnclust `out' taxi `=e(h_l)'
			addconmean `out' taxi `=e(h_l)'
			eststo
			local modelnames `modelnames' Basic
			local h = e(h_l)


			//Both [Preferred spec]
			decollinear INDdist2-INDdist17 INDpcn2-INDpcn37 if abs(taxi) < `h'
			rdrobust `out' taxi, vce(cluster pid) h(`h') covs(`r(varlist)')
			estadd local distfe = "X"
			estadd local pcfe = "X"
			estadd local metric = "1-Norm"
			estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
			addnclust `out' taxi `h'
			addconmean `out' taxi `h'
			eststo
		restore

	}
}



/************************************************************************
Votes
************************************************************************/


use "dta/crosscut_station.dta", clear

/*
Define dummy for whether the water depth---how far you have to dig to get to
water---is above median
*/
gen high_water = min_nearestdepth < 8 if !mi(min_nearestdepth)
gen low_water = min_nearestdepth >= 8 if !mi(min_nearestdepth)

//Stick to 3m around the 8m depth cutoff
keep if abs(min_nearestdepth - 8) < 3

//Keep only polling stations within 5km of a monitoring well
keep if min_nearestdist < 5



local out frelAITC

foreach sample of varlist low_water high_water {
	//Basic specification
	rdrobust frelAITC taxi if `sample'==1, vce(cluster pid)
	local h = e(h_l)
	estadd local metric = "1-Norm"
	estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
	addnclust `out' taxi `h' "`sample'==1"
	addconmean `out' taxi `h' "`sample'==1"
	eststo
	
	
	//Controlling for fixed-effects (use bandwidth of earlier specification)
	decollinear INDdist2-INDdist17 INDpcn2-INDpcn37 if `sample'==1  & abs(taxi) < `h'

	rdrobust frelAITC taxi if `sample'==1, h(`h') covs(`r(varlist)') vce(cluster pid)
	estadd local distfe = "X"
	estadd local pcfe = "X"
	estadd local metric = "1-Norm"
	estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
	addnclust `out' taxi `h' "`sample'==1"
	addconmean `out' taxi `h' "`sample'==1"
	eststo  
}






#d;
esttab using "${outpath}/results/checks/8m/main_table.tex", 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) replace nonote
sca("N Total Obs" "N_eff Obs in BW" "N_clust Clusters in BW" "conmean Control Mean" "h_l Bandwidth" "pv_rb Robust p-val"  "distfe District FEs" "pcfe Constituency FEs") 
sfmt("%9.0g" "%9.0g" "%9.0g" "%9.2f" "%9.3f") noobs
mtitles("NREGA" "NREGA" "NREGA" "NREGA""Votes" "Votes" "Votes" "Votes" ) 
	mgroups(
		"Water-Stressed" "Non-Stressed" "Water-Stressed" "Non-Stressed",
		pattern(1 0 1 0 1 0 1 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
varlabel(RD_Estimate "RD Estimate")
;
#d cr










































/************************************************************************
*************************************************************************
*************************************************************************

Appendix

*************************************************************************
*************************************************************************
************************************************************************/




/************************************************************************
*************************************************************************
Bunching
*************************************************************************
*************************************************************************/


/*
Requires DCdensity.ado to be put in your ado folder
*/
use "dta/nregs_gp.dta", clear
keep if y == 13 
keep if merge_groundwater == 3
//Histogram
hist taxi if abs(taxi) < .5, start(-.5) width(.01) xline(0, lc(black) ) color(gs10%30) ///
	graphregion(color(white)) name(hist) xlabel(-.5(.25).5, labsize(vlarge) ) ///
	xtitle("1-Norm (Taxi) Distance", size(huge))  ytitle("Density", size(huge))  


preserve
DCdensity taxi if abs(taxi) < .5, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) nograph //b(`binsize')
local est 	= string(r(theta),"%9.2f")
local se 	= string(r(se),"%9.2f")
local pval	= string(  2*( 1 - abs (   normal( r(theta) / r(se) )   ) ) , "%9.2f" )
restore

preserve
#delimit;
DCdensity taxi if abs(taxi) < .5, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //b(`binsize')
				graphopt(`" 
				xtitle("1-Norm (Taxi) Distance", size(huge))
				graphregion(color(white))
				text(4 -.4 "Log Difference:", size(vlarge) ) text(4 -.12 "`est'", size(vlarge)) 
				text(3.5 -.12 "(`se')", size(vlarge)) xlabel(-.5(.25).5, labsize(vlarge) ) name(mccrary)
				"') 
				;
								/*
				ylabel(0(.02).08, labsize(large) )
				xtitle(${runlabmaj}, size(huge) )
				ytitle("Probability Density", size(huge) ) title("`title_`eltype''", size(huge) )

				name(g`eltype')
				*/
#delimit cr

restore




graph combine hist mccrary, rows(1) graphregion(color(white) fcolor(white)) xsize(8)
graph drop hist mccrary
	
cap mkdir "results/checks/bunching"
graph export "results/checks/bunching/density_tests.pdf", replace





/************************************************************************
*************************************************************************
Balance
*************************************************************************
************************************************************************/
cap mkdir "${outpath}/results"
cap mkdir "${outpath}/results/checks"
cap mkdir "${outpath}/results/checks/balance"
cap mkdir "${outpath}/results/checks/bunching"

/************************************************************************
Station-Level Outcomes
************************************************************************/

use "dta/crosscut_station.dta", clear
rename ret_newid newid
merge 1:1 newid using "dta/prior_election_results_by_newid.dta"
/*
    Result                           # of obs.
    -----------------------------------------
    not matched                         5,726
        from master                        28  (_merge==1)
        from using                      5,698  (_merge==2)

    matched                            16,555  (_merge==3)
    -----------------------------------------
*/
drop if _m==2
rename _m merge_prior_results




eststo clear
foreach out of varlist fracAITC2011 fracAITC2009 idw_pmrb13 idw_pmrb14 {
	rdrobust `out' taxi, vce( cluster pid)
	local h = e(h_l)
	estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
	addconmean `out' taxi `h'
	addnclust `out' taxi `h'
	eststo
}

#d;
esttab , 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) noobs
sca("N Total Obs" "N_eff Obs in BW" "N_clust Clusters in BW" "h_l Bandwidth" "pv_rb Robust p-val" "conmean Mean Left of Cutoff") 
sfmt("%9.0g" "%9.0g" "%9.0g" "%9.3f") replace
mtitles("2011 AITC Vote Share" "2009 AITC Vote Share" "2013 Dry Season Drill Depth" "2014 Dry Season Drill Depth")
varlabel(RD_Estimate "RD Estimate") nonote
;
#d cr

#d;
esttab using "results/checks/balance/balance_pslev.tex", 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) noobs
sca("N Total Obs" "N_eff Obs in BW" "N_clust Clusters in BW" "h_l Bandwidth" "pv_rb Robust p-val" "conmean Mean Left of Cutoff") 
sfmt("%9.0g" "%9.0g" "%9.0g" "%9.3f") replace
mtitles("2011 AITC Vote Share" "2009 AITC Vote Share" "2013 Dry Season Drill Depth" "2014 Dry Season Drill Depth")
varlabel(RD_Estimate "RD Estimate") nonote
;
#d cr



/************************************************************************
GP-Level Outcomes
************************************************************************/

cap program drop gplevbal
program define gplevbal
	args out
	rdrobust `out' taxi 
	local h = e(h_l)
	estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
	addconmean `out' taxi `h'
	addnclust `out' taxi `h' 
	eststo
end
use "dta/nregs_gp.dta", clear

merge m:1 pid using "dta/secc_data_with_pids.dta" 
destring percent*, ignore("%") replace

keep if y==13


eststo clear
foreach out of varlist totalhouseholds percent1 percent2 percent3 percent4  {
	gplevbal `out'
}

#d;
esttab , 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) noobs
sca("N Total Obs" "N_eff Obs in BW" "h_l Bandwidth" "pv_rb Robust p-val" "conmean Mean Left of Cutoff") 
sfmt("%9.0g" "%9.0g" "%9.3f") replace
mtitles("Total HHs" "Landless Laborers" "Non-Ag Business" "Paying Income Tax" "Destitute")
	mgroups(
		"" "Percentage of HHs...",
		pattern(1 1 0 0 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
varlabel(RD_Estimate "RD Estimate") nonote
;
#d cr

#d;
esttab using "results/checks/balance/balance_gplev1.tex", 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) noobs
sca("N Total Obs" "N_eff Obs in BW" "h_l Bandwidth" "pv_rb Robust p-val" "conmean Mean Left of Cutoff") 
sfmt("%9.0g" "%9.0g" "%9.3f") replace
mtitles("Total HHs" "Landless Laborers" "Non-Ag Business" "Paying Income Tax" "Destitute")
	mgroups(
		"" "Percentage of HHs...",
		pattern(1 1 0 0 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
varlabel(RD_Estimate "RD Estimate") nonote
;
#d cr






eststo clear
foreach out of varlist percent5 percent6 percent7 percent8 percent9 percent10  {
	gplevbal `out'
}

#d;
esttab , 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) noobs
sca("N Total Obs" "N_eff Obs in BW" "h_l Bandwidth" "pv_rb Robust p-val" "conmean Mean Left of Cutoff") 
sfmt("%9.0g" "%9.0g" "%9.3f") replace
mtitles("Govt" "Public" "Private" "Below 5k Rs" "5k to 10k Rs" "Over 10k Rs")
	mgroups(
		"Percentage of HHs with Salaried Job in..." "Percentage of HHs with Monthly Income...",
		pattern(1 0 0 1 0 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
varlabel(RD_Estimate "RD Estimate") nonote
;
#d cr


#d;
esttab using "results/checks/balance/balance_gplev2.tex"
, 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) noobs
sca("N Total Obs" "N_eff Obs in BW" "h_l Bandwidth" "pv_rb Robust p-val" "conmean Mean Left of Cutoff") 
sfmt("%9.0g" "%9.0g" "%9.3f") replace
mtitles("Govt" "Public" "Private" "Below 5k Rs" "5k to 10k Rs" "Over 10k Rs")
	mgroups(
		"Percentage of HHs with Salaried Job in..." "Percentage of HHs with Monthly Income...",
		pattern(1 0 0 1 0 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
varlabel(RD_Estimate "RD Estimate") nonote
;
#d cr







eststo clear
foreach out of varlist percent24 percent25 percent26 rabi_TOTmeandaysworked  {
	gplevbal `out'
}

#d;
esttab , 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) noobs
sca("N Total Obs" "N_eff Obs in BW" "h_l Bandwidth" "pv_rb Robust p-val" "conmean Mean Left of Cutoff") 
sfmt("%9.0g" "%9.0g" "%9.3f") replace
mtitles("Unirrigated Land" "Irrigated Land" "Other Land" "2013 Avg. NREGA Days")
	mgroups(
		"Percentage of HHs Owning..." "",
		pattern(1 0 0 1 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
varlabel(RD_Estimate "RD Estimate") nonote
;
#d cr



#d;
esttab using "results/checks/balance/balance_gplev3.tex", 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) noobs
sca("N Total Obs" "N_eff Obs in BW" "h_l Bandwidth" "pv_rb Robust p-val" "conmean Mean Left of Cutoff") 
sfmt("%9.0g" "%9.0g" "%9.3f") replace
mtitles("Unirrigated Land" "Irrigated Land" "Other Land" "2013 Avg. NREGA Days")
	mgroups(
		"Percentage of HHs Owning..." "",
		pattern(1 0 0 1 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
varlabel(RD_Estimate "RD Estimate") nonote
;
#d cr









/************************************************************************
*************************************************************************
GP-Level Regression
*************************************************************************
************************************************************************/
cap mkdir "${outpath}/results"
cap mkdir "${outpath}/results/votes"
cap mkdir "${outpath}/results/votes/gplevel"

/************************************************************************
decollinear
Takes a varlist and returns a list with a non-collinear subset
************************************************************************/
cap program drop decollinear
program decollinear, rclass
	syntax varlist [if]
	
	_rmcoll `varlist' `if'
	local decoll
	foreach term in `r(varlist)' {
		if substr("`term'",1,2) != "o." {
			local decoll `decoll' `term'
		}
	}
	
	return clear
	return local varlist `decoll'
end


/************************************************************************
addnclust
Takes a variable name and calculates
************************************************************************/
cap program drop addnclust
program addnclust
	args var rv h ifstate
	if "`ifstate'" != "" {
		gen SAMPDEF = `ifstate'
	}
	else {
		gen SAMPDEF = 1
	}
	
	egen nmiss = rowmiss(`var' pid `rv')
	egen tag = tag(pid) if nmiss==0 & abs(`rv') < `h' & SAMPDEF==1
	count if tag == 1
	estadd scalar N_clust = r(N)
	drop nmiss tag SAMPDEF
end



/************************************************************************
addconmean
Takes an outcome, a variable name, and a bandwidth to calci;ate amd add
the ``control mean''
************************************************************************/
cap program drop addconmean
program addconmean
	args var rv h ifstate
	if "`ifstate'" != "" {
		gen SAMPDEF = `ifstate'
	}
	else {
		gen SAMPDEF = 1
	}
	
	egen nmiss = rowmiss(`var' pid `rv')
	sum `var' if nmiss==0 & `rv' < 0 & abs(`rv') < `h'/10 & SAMPDEF==1
	estadd scalar conmean = r(mean)
	drop nmiss SAMPDEF
end









use "dta/crosscut_gp.dta", clear

/*
Define dummy for whether the water depth---how far you have to dig to get to
water---is above median
*/
local aggvar medidw_pmr14
sum `aggvar', detail
gen low_water = `aggvar' > r(p50) if `aggvar' < . 
gen high_water = `aggvar' <= r(p50) if `aggvar' < . 
gen overall = 1

local out frelAITC
eststo clear
foreach sample of varlist overall low_water high_water {
	preserve
		keep if `sample' == 1
		//Basic specification
		rdrobust frelAITC taxi if `sample'==1
		local h = e(h_l)
		estadd local metric = "1-Norm"
		estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
		addnclust `out' taxi `h'
		addconmean `out' taxi `h'
		eststo
		
		
		//Controlling for fixed-effects (use bandwidth of earlier specification)
		decollinear INDdist2-INDdist17 INDpcn2-INDpcn37 if `sample'==1  & abs(taxi) < `h'

		rdrobust frelAITC taxi if `sample'==1, h(`h') covs(`r(varlist)') 
		estadd local distfe = "X"
		estadd local pcfe = "X"
		estadd local metric = "1-Norm"
		estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
		addnclust `out' taxi `h'
		addconmean `out' taxi `h'
		eststo
	restore
}


#d;
esttab , 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) 
sca("N_eff Obs in BW" "N_clust Clusters in BW" "h_l Bandwidth" "pv_rb Robust p-val"   "metric Metric" "distfe District FEs" "pcfe Constituency FEs") 
sfmt("%9.0g" "%9.0g" "%9.3f") replace nonote
nomtitles noobs
	mgroups(
		"Full Sample" "Water-Stressed" "Non-Stressed",
		pattern(1 0 1 0 1 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
varlabel(RD_Estimate "RD Estimate")
;
#d cr

#d;
esttab using  "${outpath}/results/votes/gplevel/app_mainresult.tex",  
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) 
sca("N Total Obs" "N_eff Obs in BW" "h_l Bandwidth" "pv_rb Robust p-val"  "distfe District FEs" "pcfe Constituency FEs") 
sfmt("%9.0g" "%9.0g" "%9.3f") replace nonote
nomtitles noobs
	mgroups(
		"Full Sample" "Water-Stressed" "Non-Stressed",
		pattern(1 0 1 0 1 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
varlabel(RD_Estimate "RD Estimate")
;
#d cr





/************************************************************************
*************************************************************************
Bandwidth Robustness
*************************************************************************
************************************************************************/
cap mkdir "${outpath}/results"
cap mkdir "${outpath}/results/checks"
cap mkdir "${outpath}/results/checks/bwfigs"

/************************************************************************
decollinear
Takes a varlist and returns a list with a non-collinear subset
************************************************************************/
cap program drop decollinear
program decollinear, rclass
	syntax varlist [if]
	
	_rmcoll `varlist' `if'
	local decoll
	foreach term in `r(varlist)' {
		if substr("`term'",1,2) != "o." {
			local decoll `decoll' `term'
		}
	}
	
	return clear
	return local varlist `decoll'
end


/************************************************************************
Vote Results
************************************************************************/

cap mkdir "estimates"
cap mkdir "estimates/bw"


//Results file
cap postclose voteest
postfile voteest h b ciu cil str10 sample str10 spec using "estimates/bw/voteest.dta", replace every(4)



use "dta/crosscut_station.dta", clear

/*
Define dummy for whether the water depth---how far you have to dig to get to
water---is above median
*/
sum idw_pmrb14, detail
gen low_water = idw_pmrb14 > r(p50) if idw_pmrb14 < . 
gen high_water = idw_pmrb14 <= r(p50) if idw_pmrb14 < . 
gen overall = 1


foreach sample of varlist overall low_water high_water {
	foreach h of numlist 0.05(.05)1 {
		//Basic specification
		rdrobust frelAITC taxi if `sample'==1, vce(cluster pid) h(`h')
		
		post voteest (`h') (`e(tau_cl)') (`e(ci_r_cl)') (`e(ci_l_cl)') ("`sample'") ("basic")

		
		//Controlling for fixed-effects (use bandwidth of earlier specification)
		decollinear INDdist2-INDdist17 INDpcn2-INDpcn37 if `sample'==1  & abs(taxi) < `h'

		rdrobust frelAITC taxi if `sample'==1, h(`h') covs(`r(varlist)') vce(cluster pid)
		post voteest (`h') (`e(tau_cl)') (`e(ci_r_cl)') (`e(ci_l_cl)') ("`sample'") ("fe")
	}

}

postclose voteest





use "estimates/bw/voteest.dta", clear

rename (ciu cil) =95
foreach v in ciu cil {
	gen `v'90 = (`v'95-b)*invnormal(.95)/invnormal(.975) + b
}


foreach sample in low_water high_water {
	foreach spec in basic fe {
		
		local title = strproper("`spec': ") + strproper(  subinstr("`sample'", "_", " ", 1) )
		
		#d;
		twoway 
			(rspike ciu95 cil95 h if spec=="`spec'" & sample=="`sample'", color(gs12)) 
			(rspike ciu90 cil90 h if spec=="`spec'" & sample=="`sample'", color(gs8)) 
			(scatter b h if spec=="`spec'" & sample=="`sample'", color(black))
			,
			graphregion(color(white)) legend(off) ylabel(-0.04(0.04)0.08)
			yline(0, lc(black)) title("`title'") name(`sample'_`spec')
			xtitle("Bandwidth")
		;
		#d cr
	}
}

graph combine low_water_basic low_water_fe high_water_basic high_water_fe, graphregion(color(white) fcolor(white)) 
graph drop low_water_basic low_water_fe high_water_basic high_water_fe

graph export "${outpath}/results/checks/bwfigs/voteest.pdf", replace

/************************************************************************
NREGA Results
************************************************************************/

cap mkdir "estimates"
cap mkdir "estimates/bw"


//Results file
cap postclose nregaest
postfile nregaest h b ciu cil str10 sample str10 spec using "estimates/bw/nregaest.dta", replace every(4)



use "dta/nregs_gp.dta", clear

keep if y > 13
local aggvar meanidw_pmrb14
sum `aggvar', detail
gen low_water = `aggvar' > r(p50) if `aggvar' < . 
gen high_water = `aggvar' <= r(p50) if `aggvar' < . 
gen overall = 1

local out rabi_TOTmeandaysworked


foreach sample of varlist overall low_water high_water {
	preserve
		keep if `sample'==1
		foreach h of numlist 0.05(.05)1 {
			//Basic specification
			rdrobust `out' taxi, vce(cluster pid) h(`h')
			
			post nregaest (`h') (`e(tau_cl)') (`e(ci_r_cl)') (`e(ci_l_cl)') ("`sample'") ("basic")

			
			//Controlling for fixed-effects (use bandwidth of earlier specification)
			decollinear INDdist2-INDdist17 INDpcn2-INDpcn37 if `sample'==1  & abs(taxi) < `h'

			rdrobust `out' taxi if `sample'==1, h(`h') covs(`r(varlist)') vce(cluster pid)
			post nregaest (`h') (`e(tau_cl)') (`e(ci_r_cl)') (`e(ci_l_cl)') ("`sample'") ("fe")
		}
		
	restore

}

postclose nregaest






use "estimates/bw/nregaest.dta", clear

rename (ciu cil) =95
foreach v in ciu cil {
	gen `v'90 = (`v'95-b)*invnormal(.95)/invnormal(.975) + b
}


foreach sample in low_water high_water {
	foreach spec in basic fe {
		
		local title = strproper("`spec': ") + strproper(  subinstr("`sample'", "_", " ", 1) )
		
		#d;
		twoway 
			(rspike ciu95 cil95 h if spec=="`spec'" & sample=="`sample'", color(gs12)) 
			(rspike ciu90 cil90 h if spec=="`spec'" & sample=="`sample'", color(gs8)) 
			(scatter b h if spec=="`spec'" & sample=="`sample'", color(black))
			,
			graphregion(color(white)) legend(off) ylabel(-1(2)5)
			yline(0, lc(black)) title("`title'") name(`sample'_`spec')
			xtitle("Bandwidth")
		;
		#d cr
	}
}

graph combine low_water_basic low_water_fe high_water_basic high_water_fe, graphregion(color(white) fcolor(white)) 
graph drop low_water_basic low_water_fe high_water_basic high_water_fe

graph export "${outpath}/results/checks/bwfigs/nregaest.pdf", replace











/************************************************************************
*************************************************************************
Other Distance Metrics
*************************************************************************
************************************************************************/
cap mkdir "${outpath}/results"
cap mkdir "${outpath}/results/checks"
cap mkdir "${outpath}/results/checks/metrics"




/************************************************************************
Votes
************************************************************************/
local labeucdist "2-Norm"
local labminswing "Inf-Norm"

foreach metric in eucdist minswing {
	use "dta/crosscut_station.dta", clear

	/*
	Define dummy for whether the water depth---how far you have to dig to get to
	water---is above median
	*/
	sum idw_pmrb14, detail
	gen low_water = idw_pmrb14 > r(p50) if idw_pmrb14 < . 
	gen high_water = idw_pmrb14 <= r(p50) if idw_pmrb14 < . 
	gen overall = 1

	eststo clear

	local out frelAITC

	foreach sample of varlist overall low_water high_water {
		//Basic specification
		rdrobust frelAITC `metric' if `sample'==1, vce(cluster pid)
		local h = e(h_l)
		estadd local metric = "`lab`metric''"
		estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
		addnclust `out' `metric' `h'
		addconmean `out' `metric' `h'
		eststo
		
		
		//Controlling for fixed-effects (use bandwidth of earlier specification)
		decollinear INDdist2-INDdist17 INDpcn2-INDpcn37 if `sample'==1  & abs(`metric') < `h'

		rdrobust frelAITC `metric' if `sample'==1, h(`h') covs(`r(varlist)') vce(cluster pid)
		estadd local distfe = "X"
		estadd local pcfe = "X"
		estadd local metric = "`lab`metric''"
		estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
		addnclust `out' `metric' `h'
		addconmean `out' `metric' `h'
		eststo
	}

	#d;
	esttab , 
	b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) 
	sca("N_eff Obs in BW" "N_clust Clusters in BW" "h_l Bandwidth" "pv_rb Robust p-val"   "metric Metric" "distfe District FEs" "pcfe Constituency FEs") 
	sfmt("%9.0g" "%9.0g" "%9.3f") replace nonote
	nomtitles noobs
	mgroups(
		"Full Sample" "Water-Stressed" "Non-Stressed",
		pattern(1 0 1 0 1 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
	varlabel(RD_Estimate "RD Estimate")
	;
	#d cr
	

	#d;
	esttab using "${outpath}/results/checks/metrics/`metric'_votes.tex", 
	b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) 
	sca("N Total Obs" "N_eff Obs in BW" "N_clust Clusters in BW" "h_l Bandwidth" "pv_rb Robust p-val" "metric Metric"  "distfe District FEs" "pcfe Constituency FEs") 
	sfmt("%9.0g" "%9.0g" "%9.0g" "%9.3f") replace nonote
	nomtitles noobs
		mgroups(
			"Full Sample" "Water-Stressed" "Non-Stressed",
			pattern(1 0 1 0 1 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
			span erepeat(\cmidrule(lr){@span})
			) 
	varlabel(RD_Estimate "RD Estimate")
	;
	#d cr


}




/************************************************************************
NREGA
************************************************************************/
local labeucdist "2-Norm"
local labminswing "Inf-Norm"

use "dta/nregs_gp.dta", clear
keep if y > 13

local aggvar meanidw_pmrb14
sum `aggvar', detail
gen low_water = `aggvar' > r(p50) if `aggvar' < . 
gen high_water = `aggvar' <= r(p50) if `aggvar' < . 
gen overall = 1

foreach metric in eucdist minswing {
	eststo clear
	local out rabi_TOTmeandaysworked
	foreach sample of varlist overall low_water high_water {
		quietly {
			
			preserve
				keep if `sample'==1
				

				//Basic specification [we'll use the bandwidth for this in all specifications]
				rdrobust `out' `metric', vce(cluster pid)
				estadd local metric = "`lab`metric''"
				estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
				addnclust `out' `metric' `=e(h_l)'
				addconmean `out' `metric' `=e(h_l)'
				eststo
				local modelnames `modelnames' Basic
				local h = e(h_l)


				//Both [Preferred spec]
				decollinear INDdist2-INDdist17 INDpcn2-INDpcn37 if abs(`metric') < `h'
				rdrobust `out' `metric', vce(cluster pid) h(`h') covs(`r(varlist)')
				estadd local distfe = "X"
				estadd local pcfe = "X"
				estadd local metric = "`lab`metric''"
				estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
				addnclust `out' `metric' `h'
				addconmean `out' `metric' `h'
				eststo
			restore

		}
	}
	#d;
	esttab, 
	b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) 
	sca("N_eff Obs in BW" "N_clust Clusters in BW" "conmean Control Mean" "h_l Bandwidth" "pv_rb Robust p-val"  "metric Metric" "distfe District FEs" "pcfe Constituency FEs") 
	sfmt("%9.0g" "%9.0g" "%9.2f" "%9.3f") mtitles("Overall" "Overall" "Water-Stressed" "Water-Stressed" "Non-Stressed"  "Non-Stressed" ) noobs
	;
	#d cr



	#d;
	esttab using "${outpath}/results/checks/metrics/`metric'_nrega.tex", 
	b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) replace nonote
	sca("N Total Obs" "N_eff Obs in BW" "N_clust Clusters in BW" "conmean Control Mean" "h_l Bandwidth" "pv_rb Robust p-val" "metric Metric"  "distfe District FEs" "pcfe Constituency FEs") 
	sfmt("%9.0g" "%9.0g" "%9.0g" "%9.2f" "%9.3f") noobs
	nomtitles
		mgroups(
			"Full Sample" "Water-Stressed" "Non-Stressed",
			pattern(1 0 1 0 1 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
			span erepeat(\cmidrule(lr){@span})
			) 
	varlabel(RD_Estimate "RD Estimate")
	;
	#d cr
}





/************************************************************************
*************************************************************************
Balance around 8m Cutoff
*************************************************************************
************************************************************************/

use "dta/nregs_gp.dta", clear
keep if y == 14
merge 1:1 pid using "dta/census2011_gp_stats.dta"

keep if min_nearestdist < 5

gen depthrec = min_nearestdepth-8

eststo clear
foreach out of varlist pc11_vd_t_hh pc11_vd_t_p pc11_vd_sc_p pc11_vd_st_p pc11_vd_rd_mdr  pc11_vd_p_sch_gov pc11_vd_p_sch_priv pc11_vd_int_caf_csc {
	rdrobust `out' depthrec
	local h = e(h_l)
	estadd scalar N_eff =  `e(N_h_r)' + `e(N_h_l)'
	sum `out' if taxi < 0 & taxi > -`h'/10
	estadd scalar conmean = r(mean)
	eststo
}


#d;
esttab , 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) 
sca("N_eff Obs in BW" "h_l Bandwidth" "pv_rb Robust p-val" "conmean Mean Left of Cutoff" /*"metric Metric" "distfe District FEs" "pcfe Constituency FEs"*/) 
sfmt("%9.0g" "%9.1f" "%9.3f" "%9.2f") replace
mtitles("Households" "Population" "SC Pop" "ST Pop" "Road" "Primary Schools" "Private Schools" "Internet Cafe")
varlabel(RD_Estimate "RD Estimate") nonote
;
#d cr



cap mkdir "${outpath}/results"
cap mkdir "${outpath}/results/checks"
cap mkdir "${outpath}/results/checks/8m"

#d;
esttab using "${outpath}/results/checks/8m/balance_at_8m.tex", 
b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01) 
sca("N_eff Obs in BW" "h_l Bandwidth" "pv_rb Robust p-val" "conmean Mean Left of Cutoff" /*"metric Metric" "distfe District FEs" "pcfe Constituency FEs"*/) 
sfmt("%9.0g" "%9.1f" "%9.3f" "%9.2f") replace
mtitles("Households" "Population" "SC Pop" "ST Pop" "Road" "Primary Schools" "Private Schools" "Internet Cafe")
varlabel(RD_Estimate "RD Estimate") nonote
;
#d cr





/************************************************************************
*************************************************************************
Is it co-partisan efficiencies?
*************************************************************************
************************************************************************/


/************************************************************************
Define programs
************************************************************************/


/************************************************************************
addnclust
Takes a variable name and calculates
************************************************************************/
cap program drop addnclust
program addnclust
	args var rv h ifstate
	if "`ifstate'" != "" {
		gen SAMPDEF = `ifstate'
	}
	else {
		gen SAMPDEF = 1
	}
	
	egen nmiss = rowmiss(`var' pid `rv')
	egen tag = tag(pid) if nmiss==0 & abs(`rv') < `h' & SAMPDEF==1
	count if tag == 1
	estadd scalar N_clust = r(N)
	drop nmiss tag SAMPDEF
end


/************************************************************************
addconmean
Takes an outcome, a variable name, and a bandwidth to calci;ate amd add
the ``control mean''
************************************************************************/
cap program drop addconmean
program addconmean
	args var rv h ifstate
	if "`ifstate'" != "" {
		gen SAMPDEF = `ifstate'
	}
	else {
		gen SAMPDEF = 1
	}
	
	egen nmiss = rowmiss(`var' pid `rv')
	sum `var' if nmiss==0 & `rv' < 0 & abs(`rv') < `h'/10 & SAMPDEF==1
	estadd scalar conmean = r(mean)
	drop nmiss SAMPDEF
end




/************************************************************************
crosscut_vilreg
Lifted from the Political Organizations paper. Does a within-GP 
difference-in-discontinuity regression on one or two "crosscut" variables
************************************************************************/
cap program drop crosscut_vilreg
program crosscut_vilreg
	syntax varlist(min=2 max=3) [if]
	marksample touse
	
	local out : word 1 of `varlist'
	local crosscut : word 2 of `varlist'
	local crosscut2 : word 3 of `varlist'
	preserve
		keep if `touse'
		quietly {
			//Basic specification [we'll use the bandwidth for this in all specifications]
			rdrobust `out' taxi, vce(cluster pid)
			local h = e(h_l)
			

			gen insamp = abs(taxi) < `h'
			gen weight = 1 - abs(taxi)/`h' if insamp
			
			if mi("`crosscut2'") {
				#d;
					reghdfe `out'  
					`crosscut'
					c.aitcmajority##c.taxi 
					c.taxi#c.(`crosscut')
					c.aitcmajority#c.(`crosscut')
					if insamp [aw=weight], cluster(pid) absorb(pid)
					;
				#d cr
			}
			else  {
				#d;
					reghdfe `out'  
					`crosscut'
					c.aitcmajority##c.taxi 
					c.taxi#c.(`crosscut')
					c.aitcmajority#c.(`crosscut')
					`crosscut2'
					c.taxi#c.(`crosscut2')
					c.aitcmajority#c.(`crosscut2')
					if insamp [aw=weight], cluster(pid) absorb(pid)
					;
				#d cr
			}
			estadd local gpfe = "X"
			estadd local metric = "1-Norm"
			estadd scalar N_eff =  e(N)
			estadd scalar h = `h'
			addconmean `out' taxi `h'
			
			if regexm("`if'", "2014") {
				estadd local year = "2014"
			}
			else {
				estadd local year = "2014-2016"
			}
			
			eststo
			drop weight insamp
		}
	restore
end



/************************************************************************
Combine datasets
************************************************************************/


/******************
Get "low water definitions"
*******************/
use "dta/nregs_gp.dta", clear

local aggvar meanidw_pmrb14
sum `aggvar', detail
gen low_water = `aggvar' > r(p50) if `aggvar' < . 
gen high_water = `aggvar' <= r(p50) if `aggvar' < . 

gen high_water_min = min_nearestdepth < 8 if !mi(min_nearestdepth) & abs(min_nearestdepth - 8) < 3
gen low_water_min = min_nearestdepth >= 8 if !mi(min_nearestdepth) & abs(min_nearestdepth - 8) < 3

gen overall = 1


keep pid low_water high_water overall low_water_min high_water_min
duplicates drop

tempfile water
save `water'





/******************
Get village-level NREGA allocations
*******************/
use "dta/village_level_nrega_and_votes.dta", clear

egen vid = group(panchayat block district village)

//A few cases where there's a missing village name
drop if mi(vid)
#d;
keep 
pid vid 
taxi euc minswing aitcmajority
rabi_meandaysworked14 rabi_meandaysworked15 rabi_meandaysworked16 
fullyear_meandaysworked14 fullyear_meandaysworked15 fullyear_meandaysworked16 
IND*
fracAITC2011 vilmean_aitc_candshare anyaitcwinner
;
#d cr

reshape long rabi_meandaysworked fullyear_meandaysworked, i(vid) j(y)


merge m:1 pid using `water'

/*


    Result                      Number of obs
    -----------------------------------------
    Not matched                             0
    Matched                            77,718  (_merge==3)
    -----------------------------------------
*/
drop if _m != 3
drop _m




gen allyears = 1
gen year2014 = y==14

egen meanaitc= rowmean(fracAITC2011 vilmean_aitc_candshare)






/************************************************************************
Run regressions
************************************************************************/

local cross meanaitc
local out rabi_meandaysworked




/******************
Define "low water" based on median
*******************/
eststo clear
crosscut_vilreg `out' `cross' if year2014 & overall == 1
crosscut_vilreg `out' `cross' if year2014 & overall == 1 & !mi(low_water)
crosscut_vilreg `out' `cross' if year2014 & low_water == 1
crosscut_vilreg `out' `cross' if year2014 & high_water == 1
#d;
	esttab, 
	b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01)  noobs
	keep(c.aitcmajority#c.*   ) 
	drop(c.aitcmajority#c.taxidist)
	varlabel(
		c.aitcmajority#c.meanaitc "Majority X AITC Support"
	) 
	sca("N Obs in BW" "N_clust Clusters" "h Bandwidth" "metric Metric" "gpfe Panchayat FEs") sfmt("%9.0g" "%9.0g" "%9.3f")
	
	mgroups(
		"Full Sample" "Subsample with Water Data",
		pattern(1 1 0 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
	mtitles("All" "All" "Water-Stressed" "Non-Stressed" ) 
	;
#d cr	

cap mkdir "${outpath}/results/checks/interpretation"
#d;
	esttab using "${outpath}/results/checks/interpretation/withingp_basic.tex", 
	b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01)  noobs replace
	keep(c.aitcmajority#c.*   ) 
	drop(c.aitcmajority#c.taxidist)
	varlabel(
		c.aitcmajority#c.meanaitc "Majority X AITC Support"
	) 
	sca("N Obs in BW" "N_clust Clusters" "h Bandwidth" "metric Metric" "gpfe Panchayat FEs") sfmt("%9.0g" "%9.0g" "%9.3f")
	
	mgroups(
		"Full Sample" "Subsample with Water Data",
		pattern(1 1 0 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
	mtitles("All" "All" "Water-Stressed" "Non-Stressed" )
	;
#d cr	


/******************
Define "low water" based on 8m cutoff
*******************/
eststo clear
crosscut_vilreg `out' `cross' if year2014 & overall == 1
crosscut_vilreg `out' `cross' if year2014 & overall == 1 & !mi(low_water_min)
crosscut_vilreg `out' `cross' if year2014 & low_water_min == 1
crosscut_vilreg `out' `cross' if year2014 & high_water_min == 1
#d;
	esttab, 
	b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01)  noobs
	keep(c.aitcmajority#c.*   ) 
	drop(c.aitcmajority#c.taxidist)
	varlabel(
		c.aitcmajority#c.meanaitc "Majority X AITC Support"
	) 
	sca("N Obs in BW" "N_clust Clusters" "h Bandwidth" "metric Metric" "gpfe Panchayat FEs") sfmt("%9.0g" "%9.0g" "%9.3f")
	
	mgroups(
		"Full Sample" "Subsample with Water Depth 8m +/- 3m",
		pattern(1 1 0 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
	mtitles("All" "All" "Water-Stressed" "Non-Stressed" ) 
	;
#d cr	

eststo clear
crosscut_vilreg `out' `cross' if year2014 & overall == 1
crosscut_vilreg `out' `cross' if year2014 & overall == 1 & !mi(low_water_min)
crosscut_vilreg `out' `cross' if year2014 & low_water_min == 1
crosscut_vilreg `out' `cross' if year2014 & high_water_min == 1
#d;
	esttab using "${outpath}/results/checks/interpretation/withingp_8mcutoff.tex", 
	b(3) se(3)  star(* 0.10 ** 0.05 *** 0.01)  noobs replace
	keep(c.aitcmajority#c.*   ) 
	drop(c.aitcmajority#c.taxidist)
	varlabel(
		c.aitcmajority#c.meanaitc "Majority X AITC Support"
	) 
	sca("N Obs in BW" "N_clust Clusters" "h Bandwidth" "metric Metric" "gpfe Panchayat FEs") sfmt("%9.0g" "%9.0g" "%9.3f")
	
	mgroups(
		"Full Sample" "Subsample with Water Depth 8m +/- 3m",
		pattern(1 1 0 0 ) prefix(\multicolumn{@span}{c}{) suffix(})   
		span erepeat(\cmidrule(lr){@span})
		) 
	mtitles("All" "All" "Water-Stressed" "Non-Stressed" )
	;
#d cr	












































/************************************************************************
*************************************************************************
*************************************************************************

Non-figure in-text result: Bootstrapped difference between estimates

*************************************************************************
*************************************************************************
************************************************************************/





local dovoteshare 	1
local donrega		1

/************************************************************************
*************************************************************************
Utility programs
*************************************************************************
************************************************************************/

/************************************************************************
decollinear
Takes a varlist and returns a list with a non-collinear subset
************************************************************************/
cap program drop decollinear
program decollinear, rclass
	syntax varlist [if]
	
	_rmcoll `varlist' `if'
	local decoll
	foreach term in `r(varlist)' {
		if substr("`term'",1,2) != "o." {
			local decoll `decoll' `term'
		}
	}
	
	return clear
	return local varlist `decoll'
end




/************************************************************************
*************************************************************************
Estimate: Dif-in-Disc for vote shares
*************************************************************************
************************************************************************/
if "`dovoteshare'" == "1" {


	/************************************************************************
	Program to run bootstrap
	************************************************************************/

	cap program drop dodif
	program dodif, eclass
		local out frelAITC

		foreach sample of varlist low_water high_water {
			//Basic specification
			rdrobust frelAITC taxi if `sample'==1, vce(cluster pidnew)	
			local h = e(h_l)
			
			//Controlling for fixed-effects (use bandwidth of earlier specification)
			decollinear INDdist2-INDdist17 INDpcn2-INDpcn37 if `sample'==1  & abs(taxi) < `h'

			rdrobust frelAITC taxi if `sample'==1, h(`h') covs(`r(varlist)') vce(cluster pidnew)
			
			local coef`sample' =  _b[RD_Estimate]
		}

		local dif = `coeflow_water' - `coefhigh_water'

		ereturn scalar dif = `dif'
	end

	/************************************************************************
	Prepare for bootstrap
	************************************************************************/

	use "dta/crosscut_station.dta", clear

	/*
	Define dummy for whether the water depth---how far you have to dig to get to
	water---is above median
	*/
	sum idw_pmrb14, detail
	gen low_water = idw_pmrb14 > r(p50) if idw_pmrb14 < . 
	gen high_water = idw_pmrb14 <= r(p50) if idw_pmrb14 < . 
	gen overall = 1

	drop if mi(idw_pmrb14,pid,taxi,frelAITC)


	egen clustorig = group(pid)
	sum clustorig

	local nclust = r(max)


	/************************************************************************
	Run bootstrap
	************************************************************************/
	cap mkdir "estimates"
	cap mkdir "estimates/bootstrap"
	local rngstate = e(rngstate)
	set seed 1689454		//From random.org
	bs e(dif), cluster(pid) idcluster(pidnew) rep(1000) saving("estimates/bootstrap/voteshare_dif.dta", replace every(10) ): dodif
	estimates save "estimates/bootstrap/voteshare_dif", replace



}








/************************************************************************
*************************************************************************
Estimate: Dif-in-Disc for NREGS
*************************************************************************
************************************************************************/
if "`donrega'" == "1" {


	/************************************************************************
	Program to run bootstrap
	************************************************************************/

	cap program drop dodifnrega
	program dodifnrega, eclass
		local out rabi_TOTmeandaysworked

		foreach sample of varlist low_water high_water {
			preserve
				keep if `sample'==1
				

				//Basic specification [we'll use the bandwidth for this in all specifications]
				rdrobust `out' taxi, vce(cluster pidnew)
				local h = e(h_l)


				//Both [Preferred spec]
				decollinear INDdist2-INDdist17 INDpcn2-INDpcn37 if abs(taxi) < `h'
				rdrobust `out' taxi, vce(cluster pidnew) h(`h') covs(`r(varlist)')
			restore
			
			local coef`sample' =  _b[RD_Estimate]
		}

		local dif = `coeflow_water' - `coefhigh_water'

		ereturn scalar dif = `dif'
	end

	/************************************************************************
	Prepare for bootstrap
	************************************************************************/

	use "dta/nregs_gp.dta", clear
	drop if y==13

	local aggvar meanidw_pmrb14
	sum `aggvar', detail
	gen low_water = `aggvar' > r(p50) if `aggvar' < . 
	gen high_water = `aggvar' <= r(p50) if `aggvar' < . 
	gen overall = 1


	drop if mi(`agvar',pid,taxi,rabi_TOTmeandaysworked)




	/************************************************************************
	Run bootstrap
	************************************************************************/
	cap mkdir "estimates"
	cap mkdir "estimates/bootstrap"
	local rngstate = e(rngstate)
	set seed 1689454		//From random.org
	bs e(dif), cluster(pid) idcluster(pidnew) rep(1000) saving("estimates/bootstrap/nrega_dif.dta", replace every(10) ): dodifnrega
	estimates save "estimates/bootstrap/nrega_dif", replace
}

